#ifndef _PEAKFACTORS_H_
#define _PEAKFACTORS_H_



struct PeakFactor
{
public:

    const char* getName()
    {
        const static char name[] = "PeakFactor"; 
        return name;
    }
    
    template
    <
        typename SignalSystem
    >
    double operator()(SignalSystem& system)
    {
        using namespace std;
       
        double sum = 0;
        double splash = 0;
        
        for (int i = 0; i < system.size() ; ++i)
        {
            vector<double> rslt(system[i].size());
            
            transform(
                system[i].begin(), system[i].end(),
                rslt.begin(),
                absolute<complex<double>,double>()
            );
        
            std::binder2nd
            <
                power<int, double>
            >
            power2(
                power<int, double>(),
                2
            );
        
            splash = max(splash, *max_element(rslt.begin(), rslt.end()));
            
            transform(
                rslt.begin(), rslt.end(),
                rslt.begin(),
                power2
            );
      
            sum =
                accumulate(
                    rslt.begin(), rslt.end(),
                    sum
                );
        }
        
        sum = sqrt(sum / (system.size() * system[0].size()));
        return splash / sum;
    }
};

struct PeakFactorN
{
    const char* getName()
    {
        const static char name[] = "PeakFactorN"; 
        return name;
    }
    
    template
    <
        typename SignalSystem
    >
    double operator()(SignalSystem& system)
    {
        return 1.0 - exp(1.0 - peakFactor(system));
    }
    
    PeakFactor peakFactor;
};

template
<
    typename Compression
>
class SpecPeakN
{
public:

    const char* getName()
    {
        const static char name[] = "SpecPeakN"; 
        return name;
    }
    
    template
    <
        typename SignalSystem
    >
    double operator()(SignalSystem& system, double compression = Compression().value)
    {
        using namespace std;
        
        if (system.empty())
            return 0;
    
        vector<double> density(system[0].size() * 2);
    
        for(int i = 0; i < system.size(); ++i)
        {
            vector<complex<double>> signal(system[i].size() * 2, 0);
            
            copy(
                system[i].begin(), system[i].end(),
                signal.begin()
            );
            
            std::vector<std::complex<double>> spec = SpecPeakN::fft(signal);
        
            transform(
                spec.begin(), spec.end(),
                density.begin(), density.begin(),
                SpecPeakN::add
            );
        }
    
        transform(
            density.begin(), density.end(),
            density.begin(),
            bind2nd(divides<double>(), system.size())
        );
    
        double power =
            accumulate(
                density.begin(),
                density.end(),
                0.0
            );
        
        auto outbandBegin = 
            density.begin()
            + min(
                (double)density.size(),
                density.size() * compression / 2.0
            );
        
        auto outbandEnd =
            density.begin()
            + min(
                (double)density.size(),
                density.size() * (1.0 - compression / 2.0)
            );
            
        double bandPower =
            accumulate(
                density.begin(), outbandBegin,
                0.0
            )
            + accumulate(
                outbandEnd, density.end(),
                0.0
            );
        
        double powerPeak =
            *max_element(
                density.begin(), density.end()
            );
        
        unsigned count =
            min(
                (double)density.size(),
                density.size() * (1.0 - compression / 2.0)
            )
            - min(
                (double)density.size(),
                density.size() * compression / 2.0
            );
    
        double avrPeak = bandPower / count;
        
        return 1.0 - exp( 1.0 - powerPeak / avrPeak);
    }

private:

    static std::vector<std::complex<double> > fft(const std::vector<std::complex<double> >& signal)
    {
        using namespace std;
        
        std::vector<complex<double> > result(signal.size());

        // Create plans
        fftw_plan plan =
            fftw_plan_dft_1d(
                min(signal.size(), result.size()),
                (fftw_complex*)(&signal[0]),
                (fftw_complex*)(&result[0]),
                FFTW_FORWARD,
                FFTW_ESTIMATE
            );

        // Forward DFT
        fftw_execute(plan);
    
        // Free memory
        fftw_destroy_plan(plan);
    
        return result;
    }

    static double add(std::complex<double> x, double y)
    {
        auto n = abs(x);
        return y + n * n;
    }
};

#endif // _PICK_FACTOR_H_
